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Abstract 



We have studied several sources of systematic uncertainty in calculating the aperture 
of the High Resolution Fly's Eye experiment (HiRes) in monocular mode, primarily 
as they affect the HiRes-II site. The energy dependent aperture is determined with 
detailed Monte Carlo simulations of the air showers and the detector response. We 
have studied the effects of changes to the input energy spectrum and composition 
used in the simulation. A realistic shape of the input spectrum is used in our anal- 
ysis in order to avoid biases in the aperture estimate due to the limited detector 
resolution. We have examined the effect of exchanging our input spectrum with a 
simple power law in the "ankle" region. Uncertainties in the input composition 
are shown to be significant for energies below 10^^ cV for data from the HiRes-II 
detector. Another source of uncertainties is the choice of the hadronic interaction 
model in the air shower generator. We compare the aperture estimate for two dif- 
ferent models: QGSJetOl and SIBYLL 2.1. We also describe the implications of 
employing an atmospheric database with hourly measurements of the aerosol com- 
ponent, instead of using an average as has been used in our previously published 
measurements of the monocular spectra. 
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1 Introduction 



The High Resolution Fly's Eye experiment consists of two air fluorescence 
detectors ( "HiRes-I" and "HiRes-II" ) located in the desert of Utah. HiRes ob- 
serves ultra-high energy cosmic rays indirectly through extensive air showers, 
i.e. cascades of secondary charged particles, which are caused by interactions 
of the primary cosmic ray particles with the earth's atmosphere. In the wake 
of the air shower, excited nitrogen molecules emit fluorescence hght in the 
ultraviolet, which is collected by mirrors and projected onto clusters of pho- 
tomultiplier tubes. Detailed descriptions of the detectors can be found in [1] 
and [2]. 

The HiRes experiment aims at measuring the arrival directions, composition, 
and flux of the most energetic cosmic rays. The two detectors allow stereo- 
scopic observation of air showers, which yields the best resolution in shower 
geometry and cosmic ray energy. An advantage of data analysis in monocular 



To whom correspondence should be addressed. E-mail: Andreas.Zech@obspm.fr 



2 



mode, i.e. separate analyses of the data from each of the two detectors, Ues 
in the higher statistics that can be achieved at the high energy end due to 
the longer lifetime of the HiRes-I detector, which started operation two years 
before HiRes-II, in 1997. Monocular analysis also allows an extension of the 
observed energy range down to energies as low as ~10^^ eV due to the larger 
elevation coverage and better time resolution of the HiRes-II detector, and also 
due to the necessity of triggering only one detector. The differential flux or 
"energy spectrum" observed in monocular mode by HiRes shows a hardening 
in the flux at around 10^^-^ eV, known as the "ankle" , and a suppression of the 
flux near 10^^-^ eV, at the expected energy of the GZK flux suppression [3,4]. 
These results have been published in [5]. 

The main systematic uncertainties that are introduced in the UHECR spec- 
trum measurement with the HiRes detectors have been reported in [6]. They 
are uncertainties in the absolute phototube calibration (±10%), the fluores- 
cence yield (±10%) and the correction for "missing energy" (±5%). The latter 
refers to the energy component that is channeled mainly into neutrinos and 
does not contribute to the ionization process. Not taking into account atmo- 
spheric effects, the uncertainty in the energy scale is ±15%, which results in a 
systematic uncertainty in the flux J of ±27%. The effect on the energy scale 
of a variation of the average vertical aerosol optical depth (VAOD) by ±1 
RMS value, from 0.04 to 0.06 and 0.02, has also been described in [6]. It is 
not larger than 9% on the average. This results in a total uncertainty in the 
energy scale of ±17%. The effect of the same VAOD variation on the aperture 
leads to an average atmospheric uncertainty in the flux J of ±15%. The total 
systematic uncertainty in the measured flux adds up to 31% for each of the 
two monocular spectrum measurements. 

In this paper, we examine additional systematic uncertainties that may affect 
the calculation of the HiRes aperture in monocular mode. Since the aperture 
of an air fluorescence detector is a function of the energy of the observed 
cosmic rays, it has to be modeled carefully with detailed Monte Carlo (MC) 
simulations. The HiRes MC simulation programs use libraries of air shower 
profiles, generated at different energies with the air shower simulation program 
CORSIKA [7] and the hadronic interaction code QGSJet [8], for a realistic 
representation of the fluctuations in the observed charged particle profiles. A 
detector response MC program simulates the light emission process along the 
shower and traces the photons through the atmosphere to the telescopes of 
the two detectors, taking into account all relevant atmospheric effects. The de- 
tector optics, electronics and trigger system are modeled in great detail using 
databases that record variable detector settings, as well as density fluctua- 
tions of aerosols in the atmosphere. After performing extensive comparisons 
between simulated events and data, which allow us to verify the quality of 
our simulations (see [9]), we estimate the detector acceptance using the ratio 
of accepted MC events {u^'^) to generated MC events (/i^^) in each energy 
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bin. To correctly simulate effects stemming from the finite energy resolution 
of the detectors and their limited elevation coverage, we use a continuous en- 
ergy spectrum and a bi-modal composition based on previous measurements 
as inputs to our simulation programs. 

The differential ffux J in each energy bin is calculated as: 



where N{Ei) is the number of observed events in the energy bin and A£' is the 
bin-width. The geometrical aperture ( solid angle) used in generating 

MC events is noted by AVL, and t is the detector live-time. Through our use of 
a continuous and realistic input energy spectrum, the finite energy resolution 
of the detectors is taken into account in the acceptance Cj = -^wu- This will 
be explained in the next chapter. In the following, we refer to the product of 
the constant Ail and the acceptance as (instantaneous) aperture. 

We will first consider the effects of varying the input energy spectrum on 
the calculated aperture and thus on the measured spectrum, in Section 2. In 
Section 3, we examine the implications of exchanging the hadronic interaction 
model in the air shower generator. For this study, we replace the QGSJet 
model, which is used in our standard spectrum measurement, with the SIBYLL 
model [10]. The effect of a variation of the assumed input composition on the 
measured spectrum is presented in Section 4. Another systematic uncertainty 
that can affect the aperture estimate of the experiment is that due to variations 
in the aerosol component of the atmosphere. For the analysis of the monocular 
spectra, we used an average atmospheric description based on measurements 
with laser systems that are installed at each detector site [11]. In Section 5, 
we re-analyze the HiRes-II monocular data with a database containing hourly 
measurements of the aerosol component of the atmosphere and compare it 
to the average description in our standard analysis. Although the systematic 
studies presented in this paper have been carried out with simulation and 
reconstruction tools of the HiRes-II analysis, their results are apphcable to 
the HiRes-I spectrum measurement as well. 



2 Input Energy Spectrum Bias 

The calculation of the cosmic ray energy spectrum from the measured energy 
distribution of events is a problem of unfolding the true spectrum of cosmic 
rays at their arrival at the earth's atmosphere from the distortions introduced 
by the detector. The energy distribution provided by the detector is a convo- 
lution of the true spectrum with the detector response, i.e. the efficiency of 



J{Ei) = N{Ei) ■ 
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the detector and its finite resolution. Following the discussion in G. Cowan's 
Statistical Data Analysis [12], the problem of unfolding can be stated in the 
following way: 

M 



Here, the true energy spectrum and the measured spectrum are divided into 
M energy bins; is the number of events in bin j of the true histogram, 
Vi the expectation value of the number of events in bin i of the measured 
histogram. Rij is the response matrix, which describes the detector response in 
each energy bin. Off-diagonal elements in are due to the limited resolution 
of the detector, which distributes a fraction of events from a certain energy 
bin over adjacent bins. 

The most straightforward way of determining the real event distribution 
from the measured values is to calculate the response matrix and apply its in- 
verse to the measured distribution. Determining the response matrix requires 
knowledge of the detector resolution and acceptance, as well as a good esti- 
mate of the true spectrum. However, as Cowan shows, even with a complete 
knowledge of Rij, this method is not applicable in most cases since it leads to 
large variances in the unfolded histogram, when the resolution is large com- 
pared to the bin-width. These variances arise due to the Poisson distribution 
of the observed data around the expectation values Vi. 

In practice, the "method of correction factors" can be applied for the decon- 
volution of the measured spectrum. This is the method used in our analysis. 
The estimator fii for the true spectrum is written as: 

fii = Cr^ ■ Hi (3) 



where are the observed data and are multiplicative correction factors 
for each energy bin. These correction factors are determined with MC simula- 
tions of both the physical model under study and the complete measurement 
process. They are just the inverse of the acceptance estimate Cj, which is given 
by the ratio of accepted over generated events in the MC in each energy bin: 



The distribution of accepted events evaluated at the reconstructed 

energies E^, whereas the distribution of generated events is given as a 
function of the true (input) energies E[. Calculation of the expectation value 
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for the estimator /tj provides an expression for the bias of the method of 
correction factors. 



The bias in the estimator E[jii\ is given by the last term of Equation 5. It 
goes to zero as the estimated acceptance, Q, approaches the true acceptance 
of the experiment, ^. The more reahstic the assumptions that go into the 
MC simulation are, the smaller the bias will be. One can estimate the bias by 
varying the model used in the simulation. 

We have calculated an estimate of the bias by varying the assumed true energy 
spectrum that is used as an input to the MC. It is useful to rewrite the term 
that describes the bias in the following way: 

C;' Cr^ -l]^,= {R- 1)/., (6) 



The bias as a fraction of the real spectrum /ij can thus be calculated from 
the ratio R of the true to the estimated acceptance. For our bias estimate, 
we assumed the true acceptance ^ to be the result of a simulation using our 
best estimate of the input energy spectrum. The estimated acceptance Q was 
calculated using a simple E~^ power law for the input energy spectrum. 

Figure 1 shows the measured energy distribution for data and a MC simulation 
assuming an E~^ input spectrum. About one third of the HiRes-II data used 
in our monocular spectrum measurement published in [5] have been included 
in this comparison. As can be seen from the distributions, and more clearly 
from the ratio plot (lower panel, data divided by MC) , this choice of the input 
spectrum is not very good. The ratio is not fiat because the assumed input 
spectrum does not have a break ("ankle"). Thus, if one normalizes data and 
MC to the same total number of events, the fraction of MC events is too small 
at low energies and too large at higher energies. We have used this MC set to 
calculate the biased acceptance estimate Cj. 

The bias has been corrected in Figure 2: instead of the E^^ spectrum, we now 
use a fit to the Fly's Eye stereo spectrum [13] to determine the shape of the 
input spectrum below the "ankle", and a linear fit to the HiRes-I spectrum 
for higher energies. The spectral index of this input spectrum is -3.01 between 
10^6-5 eV and lO^^-^ eV, -3.27 between 10^^-^ eV and 10^^-^ eV, and -2.80 above 
]^q18.7 ^]-^jg g^^(^y^ position of the "ankle" is assumed to be at 10^^'^ 

eV, corresponding to the first results of the HiRes spectrum measurements 
pubhshed in [9] . Our more recent result with higher statistics in the HiRes-II 
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Fig. 1. The top panel shows the energy distribution of HiRes-II data. The histogram 
and square points show the actual data; the triangles show the MC distribution, 
assuming an input spectrum, normalized to the same total number of events. 
The bottom panel shows the ratio of the data to MC distributions from the top plot 
with a linear fit below the "ankle" . 



data-set observes the "ankle" at 10^^'^ eV. The linear fit above the "ankle" 
is extended to the highest energies. The good agreement between data and 
MC shows that this choice of input spectrum is closer to the true spectrum /ij 
given that the MC simulates all other aspects of the experiment well, which 
was shown in [9]. This MC set is used to estimate the true acceptance ^. 

The nearly fiat ratio of the data and MC distributions in Figure 2 means that 
-wu is approximately constant if one chooses a realistic input spectrum. In 

i 

this case, according to Equation 5, the energy dependence of the expectation 
value for the true spectrum is approximately given by the input spectrum 
^f^^ . Any differences in the unfolded spectrum can be fed back into the MC 
and will improve the agreement between the energy distributions in real and 
simulated data, thus reducing the bias in the spectrum calculation with the 
updated MC simulation. 



7 




16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 20.5 




16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 20.5 

logioE(eV) 

Fig. 2. The top panel shows the energy distribution of HiRes-II data. The histogram 
and square points show the actual data; the triangles show the MC distribution, 
assuming an input spectrum based on broken power law fits to the Fly's Eye stereo 
spectrum [13] and HiRes-I spectrum [6], normalized to the same total number of 
events. The bottom panel shows the ratio of the data to MC distributions from the 
top plot with a linear fit below the "ankle". 



The bias we avoid by including the "ankle" feature in the input spectrum can 
be derived from Figure 3 which shows the ratio R of acceptances for the two 
MC simulations. A kink is visible in the ratio plot in the "ankle" region, even 
though the effect is very small. The ratio increases from ~0.97 at 10^^ eV to 
~1.07 at 10^^-^ eV and then decreases to ~0.94 at 10^^-^ eV. This bias is due 
to the limited energy resolution of the detectors which spreads event energies 
over neighboring bins. It should be noted that the same random number seeds 
have been used in both MC sets to reduce statistical fluctuations in the ra- 
tio plot. In our standard analysis, we smooth the calculated acceptance by 
replacing the acceptance histogram with an appropriate fitting function. The 
remaining statistical fluctuations are then taken into account in the spectrum 
measurement. For this study, we have not applied any smoothing procedures. 
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Fig. 3. Ratio of the acceptance calculated using an input spectrum based on the 
Fly's Eye and HiRes-I measurements, as described in the text, to the acceptance 



calculated assuming an E ^ input spectrum. 



In our analysis of the HiRes-II spectrum, we used initially an input spec- 
trum, which was soon replaced by a spectrum with the shape of the Fly's Eye 
measurement [13]. This input spectrum was used for our publication of the 
monocular spectra [6], [9]. Our measurement of the position of the "ankle" and 
the spectral index above the "ankle", as described above, was used for our up- 
dated spectrum publication [5]. Adjusting the input spectrum to feature the 
"ankle" shape has helped us to significantly improve the agreement in several 
comparison plots between data and MC events, which we use to evaluate our 
simulation programs. The effect on the acceptance is rather small, as seen in 
Figure 3. Replacing the E~'^ input spectrum with an "ankle" shape led to a 
variation in the acceptance of less than 10%. Any further adjustments of the 
exact shape of the "ankle" had negligible effects on the spectrum. 

Thus far, we have not included the observed flux suppression [5] above 10^^'* 
eV in our input energy spectrum. 
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3 Hadronic Interaction Model Uncertainty 



When simulating events, we read profiles of charged particles from a library 
of air showers and simulate the light generation, propagation, and detector 
response for different shower geometries. This "shower hbrary" contains a 
large coUection of profiles in steps of 5 g/cm^ vertical atmospheric depth at 
several fixed energies and at a zenith angle of 45°. The shower profiles were 
generated with CORSIKA for proton and iron primaries and fitted with the 
Gaisser-Hillas function [14]: 

^ exp((X^„, - X)l\) (7) 



The three fit parameters, y^max (the atmospheric slant depth at shower max- 
imum), ^max (the number of charged particles at shower maximum) and A, 
were written into the library files for each shower to characterize its profile. 
Xq was fixed at -60 g/cm^. Correlations between the mean values of the fit 
parameters and the logarithm of the shower energy are used to scale shower 
profiles from the fixed energies provided in the library to the continuous en- 
ergy spectrum required in the detector response MC, as described in [9]. Our 
analysis uses Gaisser-Hillas fits to the charged particle profiles of air showers 
to estimate the ionization energy of observed and simulated showers. The inte- 
gral over the fitted profile is multiplied by a mean ionization loss rate, derived 
from simulations with CORSIKA to be 2.19 MeV/(g cm^) [15]. Before the 
integration is carried out on MC events, the particle profile has to be adjusted 
for a fraction of 10% of the primary energy that is lost due to cuts on parti- 
cles with energies below preset thresholds in CORSIKA. We then also have to 
determine the "missing energy", which does not contribute to the ionization 
process, by comparing the estimated ionization energy of the library showers 
to their known total energy. A correction for the "missing energy" is added to 
the reconstructed energies of all simulated and real events. 

Newer CORSIKA versions provide directly information on the energy deposit 
profile of the air shower, but in this study we want to apply the same al- 
gorithms used in our published analysis of the monocular spectra. We have 
verified that our method yields results consistent with the energy deposit pro- 
files. 

The physics contributing to the electromagnetic component of the air shower 
is well understood and described in detail by the EGS code [16] within the 
CORSIKA program framework. The main uncertainty in the air shower sim- 
ulation stems from our limited knowledge of the initial hadronic interactions, 
which take place at energies by far exceeding those that can be observed in the 
laboratory. In order to get an estimate of the influence of those uncertainties 
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on the calculated apertTirc, we have generated two "shower libraries" using 
two different hadronic interaction models. CORSIKA 5.61 with QGSJetOl [8] 
has been used in our standard analysis. For this study, we generated an up- 
dated shower hbrary with CORSIKA 6.022 and QGSJetOl. The second model 
we chose was SIBYLL 2.1 [10] (with the same CORSIKA version). Differences 
between the two models can be found in their predictions of the particle mul- 
tiplicity, inelasticity and the extrapolations of the hadron-air cross-section to 
ultra-high energies. A detailed comparison is given in [17]. With regard to 
the charged particle profiles we are interested in, differences can be seen in 
the mean X^^j. values and the elongation rates ( '^'dv>g{E) ) • SIBYLL showers, 
especially in the case of proton primaries, have on average larger X^ax values 
and a slightly different elongation rate, as can be seen in Figure 8 (in the next 
chapter). Another difference between the two models is shown in Figure 4. 
Our estimates of the ionization energy fraction are roughly 2% larger when 
using SIBYLL, as compared with QGSJet. 

For the estimation of the detector aperture, we follow the same procedure 
with each of the two hadronic interaction models: the fit parameters for the 
"shower libraries" are taken from Gaisser-Hillas fits to the shower profiles. 
We found that the Gaisser-Hillas function with three parameters describes 
accurately the particle profiles for showers generated with either of the two 
hadronic interaction models. We determine the MC input composition from 
HiRes/MIA and HiRes stereo measurements of the mean X^aa; as a function 
of the cosmic ray energy. At a given energy, the mean ^max of the air shower 
distribution is correlated with the average mass of the primary cosmic ray 
flux. We assume a simple bi-modal composition of protons and iron nuclei and 
determine a proton fraction by comparing the data points with the iron and 
proton estimates given by the two models. This procedure will be described 
in more detail in the following chapter. 

Since the mean X^a^; values for pure proton showers are larger in the case 
of SIBYLL, we had to re-calculate the proton fraction that corresponds to 
the data points and adjust the input composition to contain a larger fraction 
of iron showers. The difference in the proton fractions used as input to the 
MC for the two models are shown in Figure 5. The X^^aj, distributions of 
reconstructed MC events that passed all our quality requirements are shown 
in Figure 6 for the two models. The close agreement of the distributions for 
the QGSJet and SIBYLL simulations demonstrates that we place simulated 
showers at the same distribution of atmospheric depths for either model. 

In both cases, we determine the "missing energy" from a comparison of the to- 
tal shower energy to the integral of the shower profile that has been multiplied 
by the mean ionization loss rate. Instead of applying an average correction for 
proton and iron showers, we determine the correction for the fraction of sim- 
ulated proton and iron showers that were accepted in our detector response 
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Fig. 4. Estimates of ionization energy fraction derived from shower profiles. The ratio 
of ionization energy to total energy is shown for proton (upper points, in magenta) 
and iron primaries (lower points, in black) versus the logarithm of the total energy. 
The squares are results from simulations with SIBYLL, the circles correspond to 
QGSJet. 



simulation and successfully reconstructed. 



Using the same analysis procedure for each of the two hadronic interaction 
models, we did not find any significant differences in our extensive set of com- 
parisons between distributions of data and simulated events with the two MC 
sets. Figure 7 shows the ratio of the apertures that result from simulations 
using the QGSJet and SIBYLL libraries of air showers. No smoothing algo- 
rithms have been applied to the calculated acceptances. The same random 
number seeds were used for the two MC sets to reduce statistical fluctuations. 
Both the normalization, which is consistent with 1, and the zero slope of the 
fit to this ratio show that the effect is negligible compared to the statistical 
uncertainties in our data-set. We thus find that if we apply our procedure to 
estimate the detector aperture in a consistent way, the result does not de- 
pend on the chosen hadronic interaction model. This is important since the 
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Fig. 5. Proton fractions used for the MC input composition for QGSJet and SIBYLL. 

models are continuously evolving. The latest version of the QGSJet model 
(QGSJet02) for example has been shown to generate mean Xmax values closer 
to the predictions of SIBYLL 2.1 [18]. 



4 Input Composition Uncertainty 



The fraction of air showers initiated by light and heavy (i.e. proton and iron) 
cosmic rays used in our MC simulation is determined from composition mea- 
surements by the HiRes/MIA [19] [20] and HiRes Stereo [21] experiments. Air 
fluorescence detectors like HiRes can measure the cosmic ray composition as 
a function of energy in a statistical way. The atmospheric depth, Xmax, at 
which an extensive air shower reaches its maximum size depends not only on 
the energy but also on the mass of the primary cosmic ray particle. On the 
average, heavy nuclei interact higher in the atmosphere than light nuclei of 
the same energy. Nuclei break up into fragments, each of which generates a 
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Fig. 6. The top panel shows the MC distributions of reconstructed Xmax using 
QGSJet (squares) and SIBYLL (triangles). The two distributions have been nor- 
malized to cover the same area. The bottom panel shows the ratio of the two dis- 
tributions shown in the top panel. 



sub-shower, thus distributing the initial energy over several cascades. (In a 
highly simplified picture, an iron shower is approximated by the superposi- 
tion of 56 proton showers of a factor of 1/56 smaller energies.) Unfortunately, 
statistical fluctuations between shower profiles are large and do not allow an 
event-by-event determination of the cosmic ray composition. Only the X^ax 
distribution for a given energy bin can be measured and compared to model 
predictions of protons and iron nuclei. 

Figure 8 shows the measured mean Xmax together with the pure proton and 
iron estimates from different models. We have re- interpreted the HiRes/MIA 
measurement by comparing it against the QGSJetOl model [8]. We use linear 
fits to the HiRes/MIA and HiRes Stereo points to determine an energy de- 
pendent proton fraction f{E) by comparing the fits with the simulated iron 
and proton lines of QGSJetOl. In this simple bi-modal model, we derive f{E) 
from the distance of the fitted data points to the proton and iron lines. A data 
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Fig. 7. Ratio of apertures using SIBYLL 2.1 (numerator) and QGSJetOl (denomi- 
nator) with linear fit. 

point on the proton line would have an / of 1, whereas a data point in the 
middle between the two simulation lines would have an / of 0.5. The derived 
proton fraction is 0.45 at 10^^ eV, 0.80 at 10^^-^^ eV and 1.0 at 10^° eV. The 
proton fraction we determine from the measurements and insert into our MC 
simulations depends thus on a specific interaction model. However, by using 
the same model as a reference for the input composition and for the simula- 
tion of air showers, we generate events with the measured Xmax distribution 
independently of the chosen model, as was shown in the previous section. 

Here, we investigate the effect of a change in the measured mean X^ax on our 
aperture estimate. For this study, we have used the same MC programs as 
in our standard analysis, i.e. the HiRes-II detector response simulation and 
CORSIKA 5.61 with QGSJetOl for the air shower generation. The difference 
in the estimated aperture between a MC set with only iron events and a set 
with only proton events can be seen in Figure 9. At the low energy end of 
the spectrum, the aperture for iron cosmic rays is lower because iron showers 
develop higher up in the atmosphere and are more likely to lie above the 
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Fig. 8. Xmax VS. the logarithm of the energy from HiRes/MIA (stars) and HiRes 
Stereo (triangles) measurements. The predictions for QGSJet and SIBYLL are 
shown for comparison. The diamonds are simulation points calculated by D. Heck. 
The best fit to each data set is shown. A fit to the 76% of HiRes Stereo data with 
hourly atmospheric corrections is included as well. This fit has a slightly steeper 
slope. The figure is taken from [21]. 



HiRes-II elevation coverage (3° to 31°) than proton showers. This leads to 
larger differences between the two apertures at lower energies. For energies 
above ~ 10^^ eV, where showers are on average farther away from the detector, 
no significant difference is seen in the aperture for iron and proton showers. 

Using the pure iron and pure proton apertures, we have calculated the effect of 

a change in the assumed proton fraction on the aperture estimate. The proton 
fraction f{E) is defined as the ratio of generated proton showers over the sum 
of generated proton and iron showers in the MC: f{E) = ij,p{E) / {iip{E) + 
IJ,fe{E)). The acceptance for a MC set with a proton fraction / in a given 
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Fig. 9. Apertures for pure iron and pure proton MC sets. 



energy bin is: 

Up + Vfe 



a/ 



ap{R + f-{l-R)) 



Here Up and uje are the accepted, i.e. triggered and well reconstructed, proton 
and iron events, respectively; af and are the acceptances for a MC set with 
proton fraction / and 1, respectively. R is given by the ratio of the acceptances 
for pure iron and pure proton MC sets (-jr/jT^^ ■ This ratio can be determined 
directly from the two curves shown in Figure 9, since the apertures are just 
the acceptances multiplied by a constant factor, the geometrical aperture AQ. 
With R known. Equation 8 yields the acceptance a/ for a given proton fraction 
/ in a given energy bin. 

It can be seen from Figures 9 and 8 that systematic uncertainties in the aper- 
ture due to uncertainties in the proton fraction are only of concern at the low 
energies covered by the HiRes/MIA measurement. We have calculated the sys- 
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tematic uncertainty in the proton fraction f{E) from the relevant uncertainties 
in energy and Xj^ax quoted in the HiRes/MIA PRL paper [20]. 

Sources for energy uncertainties in HiRes/MIA are the detector cahbration 
(<5% uncertainty in energy) and the aerosol component of the atmosphere 
(<10% uncertainty in energy). A 10% uncertainty in the fluorescence yield 
is common to both HiRcs and HiRcs/MIA, and is therefore not included in 
our calculation. Since both experiments use the same assumptions on the 
fluorescence yield, a potential error in this parameter would induce the same 
bias in the reconstructed energies of HiRes/MIA and HiRes. It would thus 
not change the shape of the aperture. Given the measured elongation rate of 
93 g/cm^ [20], the uncertainties in energy from calibration and atmosphere, 
added in quadrature, contribute <4.4 g/cm^ to the uncertainty in X^^x- 

The quoted uncertainty in Xj^ax of roughly 25 g/cm^ due to the calculation 
of the Cherenkov fraction is also common to the two experiments and is thus 
not relevant for our calculation. Since the same assumptions on the Cherenkov 
light beam are made in the HiRes and HiRes/MIA analysis, a potential bias 
in the HiRes/MIA reconstruction would be corrected in the HiRes detector 
simulation before the calculation of the acceptance. In other words, the HiRes 
MC simulation positions showers on the average at the same height where they 
were seen by HiRes/MIA. A recent study of the fluctuations of the molecular 
density proflle using radio sonde data shows a small discrepancy with the 
standard model used in both HiRes/MIA and HiRes [22]. This introduces an 
additional uncertainty in X^ax of <10 g/cm^. 

Since the separation between the proton and iron lines in the QGSJetOl model 
is ~100 g/cm^, the uncertainties in X„iax of 4.4 g/cm^ from the energy mea- 
surement and of 10 g/cm^ from the molecular density fluctuations translate 
to ~4.4 % and ~10 % uncertainty in the proton fraction /(-£'), respectively. 
Finally, one has to add a ~3% uncertainty coming from the linear flt to the 
HiRes/MIA data that is used to parameterize the proton fraction in the sim- 
ulation programs. Those uncertainties in the HiRes/MIA measurement of the 
mean X^ax that translate into uncertainties in the acceptance add then up to 
~11% of the proton fraction. 

With the help of equation 8, we have calculated the variation in the aperture 
af{E) for a variation in f{E) of ±11%. A change in the aperture translates 
directly into a change of the measured spectrum. The uncertainties in the 
spectrum from a ±11% variation in the proton fraction are shown in Figure 10 
as thick error bars. At the low energy end of the HiRes-II spectrum, the 
systematic uncertainties from the input composition are comparable to the 
statistical uncertainties in the spectrum. In the absence of a more precise 
composition measurement in the HiRes/MIA energy range, it will thus be 
difficult to observe the feature of the "second knee" [23], [24], [25], [26], even 
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Fig. 10. HiRes-11 energy spectrum with systematic uncertainties (thick error bars) 
corresponding to a ±11% change in the proton fraction of the MC. Data included 
in the spectrum were recorded between December 1999 and September 2001. 



5 Atmospheric Uncertainties 



In the uhraviolet (UV) fluorescence range observed by the HiRes telescopes 
(310-400 nm), light attenuation comes from mainly two sources: Rayleigh scat- 
tering on air molecules, and absorption and scattering on aerosols. While 
Rayleigh scattering is well understood and does not vary much over time, 
attenuation by aerosols can change over short time ranges and has to be mon- 
itored during the data taking process. We use a system of steerable lasers, one 
at each of the two detector sites, to measure light attenuation by aerosols on an 
hourly basis. The UV laser at each site fires shots in a regular pattern of vary- 
ing geometries, which are observed from the other site. The vertical aerosol 
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optical depth (VAOD) can be measured from the detected hght of vertical 
shots. Horizontal attenuation length (HAL) and scattering phase function due 
to aerosols can be measured from the light scattered into the telescopes under 
different angles from horizontal shots. The wavelength of the laser is 355 nm, 
close to the 357 nm fluorescence line. We account for the wavelength depen- 
dence of aerosol scattering in our simulation and reconstruction programs. 

For the monocular spectra published in [5] and [6], we used a measurement 
of the average VAOD and HAL in our analysis. This was necessary since the 
steerable laser system became fully operational only two years after the HiRes- 
I detector had started taking data. For consistency in the analysis of the two 
monocular spectra, we have thus used a single average measurement, while 
applying strict cuts on the selection of clear nights that were included in the 
spectrum. Here, we repeat the analysis for HiRes-II using a database with 
hourly entries of the measured VAOD and HAL instead of the average values. 
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Fig. 11. Hourly VAOD and inverse HAL measurements for the selection of clear 
nights in the HiRes-II analysis (December 1999 to September 2001). 



The inverse HAL and VAOD distributions for the selected clear nights that 
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went into this analysis are shown in Figure 11. In the atmospheric database 
generated for the HiRes-II analysis, entries are available for 82 out of the 122 
nights that were selected for the spectrum calculation. When measurements 
of both parameters were not available, a seasonal average was assumed. The 
averages found in this analysis represent a slightly clearer atmosphere than 
the values used in the standard monocular analyses (<HAL> = 25 km and 
<VAOD> = 0.04 were used for our published spectra). This difference is due 
to the data normalization method applied here, which was not used in the 
original analysis (see [11] and [27]). The averages determined here (<HAL> 
= 27 km and <VAOD> = 0.035) are nevertheless well within the quoted 
uncertainties of the averages used in the monocular analyses. 

In order to study the effect of variations in the aerosol component of the atmo- 
sphere on the reconstructed energies, we have analyzed the HiRes-II data from 
December 1999 to September 2001 using the atmospheric database. Since all 
events were reconstructed both with the atmospheric database and with the 
average atmospheric values, the ratio of the energy estimates can be calculated 
for each event. The distribution of those ratios is shown in Figure 12 as a func- 
tion of the energy reconstructed using the atmospheric database. A Gaussian 
fit has been applied to the distribution in each energy bin. The points repre- 
sent the Gaussian means, the error bars the standard deviations. The energies 
reconstructed with database are on the average 4% smaller. This is due to the 
slightly clearer atmosphere determined with the improved analysis method for 
the VAOD values that went into the database. There is no significant energy 
dependence in the ratio of the energies. 

We then examined the effect of the small shift in reconstructed energies on the 
distribution of events over energy bins N{Ei), which goes into the spectrum 
calculation. The histograms for the two energy reconstructions, using the aver- 
age atmosphere and the atmospheric database, can be seen in Figure 13. The 
4% shift in energy is too small to cause a significant effect in the event distri- 
bution given our bin-size, which is adapted to the data statistics.lt should be 
noted that the difference between the two distributions is here even aggravated 
by the fact that two slightly different versions of our reconstruction software 
were used. The histogram for the average atmosphere is the exactly same as 
in the calculation of the published HiRes-II spectrum, which permits a direct 
comparison of this figure with Figure 15, wheres for the reconstruction with 
database we had to use a slightly updated version of our analysis software. 
(For Figure 12, only the updated version was used.) 

The effect of atmospheric variations on the energy resolution can be studied 
with simulated events. We have generated a MC set with about four times 
data statistics using the atmospheric database. The MC events have been 
reconstructed in two different ways: first with the seasonally averaged atmo- 
spheric values and then with the database for nights when atmospheric data 
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Fig. 12. Difference in reconstructed energies due to atmospheric variations. The en- 
ergy of each data event has been reconstructed using both the atmospheric database 
and an average atmosphere. Shown is the ratio of the two reconstructed energies 
(average over database) versus the energy reconstructed using the atmospheric 
database. The sizes of the boxes are proportional to the number of entries. The 
markers represent the mean and width of a Gaussian fit to the distribution in each 
energy bin. 



were available. A comparison of the resolution estimates is shown in Figure 14. 
It should be noted that the two plots use a logarithmic scale, hence the tails 
in the distributions are very small. There is no significant difference between 
the tails of the two distributions. Only the width differs by a small amount. 
Reconstructing the MC events with seasonally averaged atmospheric values in- 
stead of using the atmospheric database widens the resolution by 0.9% (both 
in RMS and a of the Gaussian fit). 

Finally, we have analyzed the effect of using the atmospheric database rather 
than the measured average on the energy spectrum. We have calculated the 
acceptance from a MC set that was generated and reconstructed with the 
database. As in our standard analysis, we have applied a smoothing proce- 
dure to minimize statistical fluctuations in the acceptance. The remaining 
statistical uncertainties of the smoothed acceptance are taken into account in 
the measurement of the spectrum. The HiRes-II data were also reconstructed 
with use of the atmospheric database. In this way, the hourly measurements of 
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Fig. 13. The top panel shows the energy distribution of HiRes-II data analyzed 
using the atmospheric database (triangles) and analyzed using average atmospheric 
parameters (squares). The bottom panel shows the ratio of the two distributions 
with the average atmosphere case as the numerator and the atmospheric database 
case as the denominator. 



atmospheric variations were included in every step of the analysis. The energy 
spectrum resulting from this analysis is compared to the published spectrum, 
which uses the nominal averages of VAOD (0.04) and HAL (25 km), in Fig- 
ure 15. The result for JE^ does not vary by more than ±15% at any energy, 
except for the first and the last two bins, where statistics in the data are lim- 
ited. The difference between the two spectra in the last two bins is due to a 
single event that has shifted up in energy to the last bin when reconstructed 
using the atmospheric database, as can be seen from Figure 13. 

The focus of the study presented in this chapter was on the difference between 
reconstruction results with an atmospheric database and with an average at- 
mosphere. We have also examined the systematic uncertainty in the determi- 
nation of the atmospheric parameters that describe the aerosol distribution 
(HAL and VAOD). We have compared the atmospheric database used in this 
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Fig. 14. The energy resolution is shown as the natural logarithm of the ratio of 
reconstructed over true (input) energy. The top panel shows the resolution in a MC 
set generated using the atmospheric database and reconstructed with seasonally 
averaged atmospheric parameters. For the ratio of reconstructed over true energy 
(not its logarithm) one obtains a width of 21.9% (RMS) and a o" of the Gaussian 
fit of 18.9%. The bottom panel shows the resolution in the same MC set when 
analyzed using the database. The RMS width and a of the ratio are 21.0% and 
18.0%, respectively. 



study with an independent result from an analysis of laser shots with different 
geometries, which were reconstructed with independently developed software. 
The reconstructed energies one obtains with this independent database are on 
the average 5% lower than the values presented here. The average difference 
remains smaller than 7% at the highest energies. 
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Fig. 15. The HiRes-II spectrum from our standard analysis is shown in circles. The 
spectrum resulting from an analysis that was using the atmospheric database in 
the determination of the aperture and in the reconstruction of the data is shown in 
triangles. 



6 Conclusions 



None of the sources of possible systematic uncertainties we have studied here 
contribute significantly to our published estimate of the systematic uncer- 
tainty. The bias introduced by using an power law instead of a more 
realistic spectral shape is not very significant in the "ankle" region. Our calcu- 
lated aperture is sensitive to the assumed input composition for energies below 
^10^*^ eV for HiRes-II. By using a measured composition as an input to our 
simulation programs, our analysis does not depend on the assumed hadronic 
interaction model. For the 17 month period tested here, the description of 
the aerosol density using an hourly database does not cause any significant 
differences in the spectrum, when compared with an average atmosphere. We 
also found no significant changes in the reconstructed energies for the time 
period under study. 



25 



Acknowledgements 



This work is supported by US NSF grants PHY-9321949, PHY-9322298, PHY- 
9904048, PHY-9974537, PHY-0098826, PHY-0140688, PHY-0245428, PHY- 
0305516, PHY-030098, and by the DOE grant FG03-92ER40732. We gratefully 
acknowledge the contributions from the technical staffs of our home institu- 
tions. The cooperation of Colonels E. Fischer and G. Harter, the US Army, 
and the Dugway Proving Ground staff is greatly appreciated. 



References 



[1] T. Abu-Zayyad et al, Proceedings of the 26th International Cosmic Ray 
Conference 5, 349 (1999) 

[2] J. Boyer, B. Knapp, E. Mannel and M. Seman, Nucl. Inst, and Meth. A 482, 
457 (2002) 

[3] K. Greisen, Phys. Rev. Lett. 16, 748 (1966). 

[4] G.T. Zatsepin and V.A. K'uzmin, Pis'ma Zh. Eksp. Teor. Fiz. 4, 114 (166) 
[JETP Lett. 4, 78 (1966)]. 



[5] R. U. Abbasi et al, Phys. Lett. B 619, 271 (2005), astro-ph/0501317 



[6] R. U. Abbasi et al, Phys. Rev. Lett. 92, 151101 (2004), [ asto>ph/02082 43 

[7] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz and T. Thouw, FZKA 6019 
(1998), Forschungszentrum Karlsruhe, 



http: / / www-ik3.fzk.de/~heck/corsika/physics_description/corsika_phys.html 



[8] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, Nucl. Phys. B (Proc. 
Suppl.) 52B, 17 (1997) 

[9] R. U. Abbasi et al, Astropart. Phys. 23, 157 (2005), [astro^h/020830T] 

[10] R. S. Fletcher, T. K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D, 50, 5710 
(1994) 

[11] R. U. Abbasi et al, Astropart. Phys. 25, 74 (2006) 

[12] G. Cowan, Statistical Data Analysis, Oxford Science Publications, 1998 

[13] D. J. Bird et al, Phys. Rev. Lett. 71, 3401 (1993) 

[14] T. K. Gaisser and A. M. Hillas, Proceedings of the 15th International Cosmic 
Ray Conference, 8, 353 (1977) 



[15] C. Song et al, Astropart. Phys. 14, 7 (2000), astro-ph/9910195 



[16] W. R. Nelson, H. Hirayama and D. W. O. Rogers, SLAC Report, 265 (1985) 



26 



[17] D. Heck, M. Risse and J. Knapp, Nucl. Phys. B (Proc. Suppl.) 122, 364 (2003) 

[18] S. S. Ostapchenko, presentation at the VIHKOS CORSIKA school 2005 in 
Lauterbad, Germany 

[19] T. Abu-Zayyad et al, Astrophys. J. 557, 686 (2001), |as tro-ph/0010652 

[20] T. Abu-Zayyad et al, Phys. Rev. Lett. 84, 4276 (2000) 

[21] R. U. Abbasi et al, Astrophys. J. 622, (2005) 910-926 

[22] Y.Fedorova (HiRes Collaboration), Proc. of the 29th Int. Cosmic Ray Conf. 
(2005) 

[23] D. J . Bird et al, Ap. J. 424, 491 (1994) 

[24] T. Abu-Zayyad et al, P.R.L. 84, 4276 (2000) 

[25] N. Nagano et al, J. Phys. G 18, 423 (1992) 

[26] M. I. Pravdin et al, Proc. 26th ICRC, Salt Lake City (1999). 

[27] R. U. Abbasi et al, Astropart. Phys. 25, 93 (2006) 



27 



